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Abstract — 

The  finite  element  (FE)  formulation  for  waveguide  discontinuity 
analysis  is  reviewed  and  extended  to  multiple,  arbitrarily-oriented 
ports.  Several  higher-order  vector  elements  — specifically  hierarchal 
linear  tangential/quadratic  normal  (LT/QN)  — are  compared,  and  the 
extensions  required  to  incorporate  LT/QN  elements  in  the  formulation 
are  presented.  The  improved  accuracy  afforded  by  LT/QN  elements 
compared  to  constant  tangential/linear  normal  (CT/LN)  elements  is 
investigated  by  considering  energy  conservation  in  an  empty  waveg- 
uide. Results  obtained  using  both  CT/LN  and  LT/QN  elements  are  also 
shown  for  a problem  of  engineering  interest:  an  E-plane  bend.  Results 
for  the  LT/QN  elements  compare  especially  well  to  approximate  ana- 
lytical results  using  quite  coarse  meshes.  The  paper  concludes  with  a 
discussion  of  the  use  of  iterative  solvers  and  possible  convergence  prob- 
lems encountered  when  using  higher-order  elements. 

Keywords — Finite  element  method;  higher-order  vector  elements; 
waveguide  discontinuities. 

I.  INTRODUCTION 

The  analysis  of  waveguide  discontinuities  has  been  a 
canonical  problem  for  analytical,  approximate,  and  now  nu- 
merical approaches  since  the  pioneering  work  of  Marcuvitz 
and  colleagues  during  the  Second  World  War,  now  some 
sixty  years  back.  Using  variational  formulations,  and  quasi- 
static approximations  of  the  fields,  Marcuvitz  et  al.  were 
able  to  analyze  an  extraordinary  variety  of  problems,  docu- 
mented in  the  classic  text  [I].  Subsequently,  mode-matching 
methods  were  introduced  for  the  analysis  of  “stepped”  dis- 
continuities — i.e.  structures  where  the  waveguide  modes 
could  be  computed  in  a step-wise  fashion,  and  matched  at 
two-dimensional  planes.  However,  for  general,  arbitrary  dis- 
continuities, and  of  course  those  involving  non-metallic  dis- 
continuities such  as  dielectrics,  differential  equation  based 
methods  such  as  the  finite  element  method  (FEM)  and  finite 
difference  time  domain  (FDTD)  method  are  now  the  meth- 
ods of  choice. 

Although  an  obvious  application  of  the  FEM,  disconti- 
nuities in  rectangular  waveguide  have  not  been  widely  ad- 
dressed in  the  literature,  in  particular  using  higher-order  el- 
ements. Ise  [2]  used  “brick”  elements  of  “first”  order  (con- 
stant tangential  / linear  normal  — CT/LN)  to  analyze  both  a 
dielectric  post  and  a concentric  step  discontinuity  in  rectan- 
gular waveguide  ; Jin  presented  a detailed  formulation  in  [3, 


Chapter  8],  also  using  CT/LN  elements;  Webb’s  review  pa- 
per discussed  a number  of  related  issues  but  did  not  address 
higher-order  elements  [4];  and  Pekel  and  Lee  addressed  the- 
oretical aspects  of  mesh  refinement  using  an  empty  piece  of 
waveguide,  but  again  did  not  explicitly  discuss  higher-order 
elements  [5].  Scott  addressed  rotationally  symmetric  waveg- 
uide and  obtained  very  good  results  using  special-purpose 
higher-order  elements  [6]. 

The  contributions  of  this  paper  are  the  following.  Firstly, 
Jin’s  formulation  is  extended  to  arbitrarily  oriented  waveg- 
uides (Jin’s  formulation  assumes  z orientation),  with  multi- 
ple ports  (Jin  assumes  two  port  devices).  Secondly,  the  avail- 
able higher-order  vector  elements  are  reviewed,  and  some 
unifying  themes  underlying  these  are  identified.  Thirdly,  the 
necessary  extensions  to  include  higher  order  vector  elements 
in  the  formulation  are  outlined.  Fourthly,  the  accuracy  ob- 
tained vs.  element  size  and  number  of  degrees  of  freedom 
for  CT/LN  and  linear  tangential  / quadratic  normal  (LT/QN) 
elements  is  investigated  by  monitoring  energy  conservation 
in  a piece  of  empty  guide.  Finally,  the  extended  formula- 
tion and  implementation  is  validated  — and  the  far  greater 
accuracy  obtainable  with  the  LT/QN  elements  demonstrated 
again  — by  analyzing  a realistic  waveguide  problem  in  X- 
band  waveguide,  namely  an  E-plane  bend.  Results  for  this 
are  compared  with  Marcuvitz ’s. 

Some  aspects  of  this  paper  were  originally  presented  in 
[7].  However,  the  formulation  presented  therein  is  Jin’s, 
and  does  not  incorporate  the  new  extensions  to  be  presented 
here,  which  are  required  to  analyze  general  waveguide  struc- 
tures (such  the  E-plane  bend  analyzed  here).  Furthermore, 
the  discussion  of  higher-order  elements  has  been  extensively 
revised,  to  highlight  the  connection  between  different  pub- 
lished elements.  Finally,  some  problems  regarding  conver- 
gence of  the  iterative  solvers  which  have  emerged  subse- 
quent to  [7]  are  discussed. 

II.  The  waveguide  formulation 

The  formulation  is  a straightforward  extension  of  Jin’s  ap- 
proach [3].  His  formulation  addressed  two-port,  single  mode 
analysis,  with  the  waveguide  oriented  in  the  ^-direction. 
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Here,  general  waveguide  orientation(s)  will  be  considered. 
The  formulation  assumes  hollow,  rectangular  guide  at  the 
ports  (although  the  extension  to  homogeneously  filled  guide 
is  straightforward).  The  TEi0  mode  is  assumed  in  the  fol- 
lowing. In  between  the  ports,  in  the  region  to  be  discretized 
using  finite  elements,  the  waveguide  may  contain  linear,  in- 
homogeneous, lossy,  dielectric  and/or  magnetic  material(s); 
and/or  conductors  (for  instance,  posts  or  irises);  and  may 
change  orientation  (eg  E-plane  bends)  or  dimension  (eg  E- 
and/or  H-plane  steps).  The  formulation  to  be  used  does, 
however,  assume  isotropic  media.  The  generalization  of 
the  analysis  to  multiple  ports,  the  inclusion  of  higher-order 
modes,  and  the  extension  to  more  general  waveguide,  will 
be  outlined  subsequently. 

A.  Formulation  overview 

The  key  part  of  the  formulation  is  to  write  the  electric  field 
at  port  1 (5 1)  as  the  sum  of  the  known  incident  and  unknown 
reflected  fields  in  terms  of  the  C)  coordinate  system 
local  to  the  port,  with  £ in  the  direction  of  propagation,  and 
set  to  zero  at  each  port,  as  follows: 

E(S,v,0  = £inc(£,»7,0 + £ref  (£,*?,  0 

= (E0eio(t,v)e~:lkc'‘0<  + 

RE0e1Q^,v)e+jk^%=0  (1) 

eio(£,r?)  is  the  relevant  waveguide  eigenmode  (the  TEi0 
eigenmode  here)  and  k^l0  is  the  modal  propagation  constant. 
Note  that  it  is  necessary  to  retain  the  term,  even 

though  the  field  is  evaluated  at  £ = 0,  since  the  boundary 
condition  to  be  discussed  involves  the  derivative  of  the  field, 
which  must  be  evaluated  before  setting  £ = 0. 

The  next  key  element  of  the  formulation  is  to  convert 
eqn.  (1)  to  a boundary  condition  of  the  third  type  involving 
both  the  field  and  its  normal  derivative;  the  detail  is  given  in 
[3,  §8.5]: 

n x ( V x E)  + 7n  x (n  x E)  = Uinc  (2) 

with 

7 = U™  = -2jk<l0Einc  (3) 

It  should  be  noted  that,  in  obtaining  eqn.  (2),  the  transverse- 
only  nature  of  the  TE  field  is  exploited.  TM  modes  contain 
axial  E field  components,  and  the  boundary  condition  cannot 
thus  be  written  for  an  E field  solver.  TM  mode  analysis 
could  be  undertaken  by  using  an  H field  solver. 

The  same  is  repeated  at  port  2,  but  at  that  port,  there  is 
only  an  unknown  transmitted  field: 


E(S,1 1,0  = &I&ns(Cv, 0 

= (4) 


Similar  comments  apply  as  regards  the  e ik<io£  term.  The 
boundary  condition  at  port  two  is 

h x ( V x E)  + 7n  x (n  x E)  = 0 (5) 

In  Jin’s  original  formulation,  eqn.  (4)  was  written  as 


E(x,y,z)  = Eu*ns(x,y,z) 

= TE0e10(x,y)e-^oz  (6) 

In  this  approach,  z — Z\zX  port  \9z  = z?  at  port  2,  thus  the 
phase  in  Jin’s  formulation  was  referenced  to  each  port.  In  the 
present  formulation,  the  transmission  coefficient  T incorpo- 
rates the  “insertion”  phase  — i.e.  for  a section  of  empty 
guide  length  f T will  have  phase  angle  - kZl0£ , whereas  in 
Jin’s  formulation,  the  phase  was  0.  The  present  formula- 
tion produces  the  same  phase  that  would  be  measured  us- 
ing a vector  network  analyzer,  with  reference  planes  cali- 
brated at  the  ports.  (Jin’s  approach  worked  well  for  straight 
waveguide,  but  is  inappropriate  for  “bent”  guides  or  multiple 
ports.) 

The  equivalent  variational  functional  (assuming  isotropic 
but  possibly  lossy  materials),  subject  to  these  boundary  con- 
ditions on  the  ports  and  Etan  = 0 on  the  perfectly  conduct- 
ing walls,  is  well  known: 


F(E) 


J j J [—(V  x E)  • (V  x E)  - kl<irE-  i?J  dV 

//j 

J J \l(n  x E)  ■ (n  x £)]  dS 


* 7-(n  x E)  ■ (n  x E)  + E • t?incj  dS 


The  FE  discretization  of  this  functional  is  discussed  in  Sec- 
tion IV. 


B.  Computation  of  the  S-parameters 

The  above  formulation  produces  R and  T for  port  1 (Sn 
and  S2i)-  It  must  be  repeated  with  an  incident  field  at  port  2 
to  obtain  S 12  and  522-  Only  the  excitation  vector  changes,  so 
this  is  simply  a question  of  repeating  the  matrix  solve.  For 
multiple  ports,  the  extension  is  obvious:  T is  computed  at 
each  port,  producing  one  column  of  the  5 matrix.  The  exci- 
tation is  then  repeated  at  each  port  to  produce  other  columns. 
Although  it  will  not  be  shown  in  this  paper,  the  formulation 
has  also  been  verified  successfully  by  the  author  for  a four- 
port  device. 

The  S-parameters  may  be  computed  directly  from  the 
fields  on  the  ports.  A more  accurate  approach  uses  the  or- 
thogonality of  the  modes  to  integrate  the  fields  computed 
over  each  port  [3,  §8.5];  as  an  example  for  the  two-port  ex- 
ample, the  transmission  coefficient  is  given  by: 

T = ^0  I Is  ^ ^ ° ' e’*lo(^ T?)d5  (8) 
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As  before,  eio(£,  rj)  is  the  relevant  waveguide  eigenmode;  a 
and  b are  the  waveguide  dimensions. 

C.  The  waveguide  formulation:  another  perspective 

The  formulation  can  be  viewed  as  a finite  ele- 
ment/boundary integral  (FE/BI)  formulation,  using  the 
waveguide  Green’s  function  for  “exact”  mesh  termination. 
(For  radiation  or  scattering  problems,  FE/BI  formulations 
use  the  free  space,  or  sometimes  the  half-space,  Green’s 
function,  eg.  [3,  §9.3]).  The  current  dominant-mode-only 
analysis  uses  only  the  first  in  the  infinite  series  of  modes 
comprising  the  waveguide  Green’s  function.  It  is  accurate 
provided  that  the  ports  are  sufficiently  far  removed  from  the 
discontinuities  (assuming,  of  course,  that  only  the  dominant 
mode  is  above  cut-off).  Higher  order  modes  are  easily  in- 
cluded in  the  formulation;  this  does  require  re-computing 
both  the  LHS  matrix  and  RHS  vector,  since  the  former  has 
one  term  dependent  on  the  propagation  constant,  and  the 
latter  is  obviously  dependent  on  the  incident  mode  shape. 
The  formulation  presently  assumes  hollow  waveguide  at  the 
ports;  i.e.  only  TE  (and  TM  modes,  if  an  H field  solver  is 
also  implemented)  are  included.  More  exotic  modes,  or  nu- 
merically determined  ones,  could  also  be  incorporated  into 
the  formulation.  This  is  discussed  further  in  Section  VI. 

III.  Higher  ORDER  VECTOR  ELEMENTS 

A.  Vector  elements  — a brief  review 

Edge  elements  (also  known  as  Nedelec  elements;  Whit- 
ney elements/forms;  CT/LN  elements;  iTo(curl)  elements) 
were  introduced  during  the  1980’s.  Nedelec  [8]  provided  the 
mathematical  framework  for  mixed  order  finite  elements  of 
various  order.  However,  the  polynomial  spaces  from  which 
the  basis  functions  were  to  be  chosen  were  defined  by  him 
in  terms  of  Cartesian  coordinates,  which  is  not  the  form  vec- 
tor elements  are  generally  given  in  now.  Cendes,  Bossavit, 
Webb  and  others  introduced  vector  elements  to  EM  FEA 
analysis  during  the  late  1980’s  and  1990’s.  (See  [9]  for  the 
original  references).  The  element  shape  function  was  then 
presented  in  terms  of  simplex  coordinates  [10],  as  what  is 
now  recognized  as  a Whitney  form,  dating  back  to  much  ear- 
lier work  by  Whitney: 

Wij  = XiVXj  ~ XjVXi  (9) 

This  element  has  the  well-known  properties  of  constant 
tangential/linear  normal  field  (CT/LN)  approximation  along 
edges  (hence,  of  mixed  order).  Since  the  approximation  is 
constant  in  the  direction  tangential  to  the  edge  connecting 
nodes  i and  j,  and  perpendicular  to  all  the  other  edges  (two, 
for  triangles,  or  five,  for  tetrahedrons),  the  degrees  of  free- 
dom, defined  by  Nedelec  as  the  line  integrals  of  the  finite 
element  approximation  along  the  respective  edges,  are  sim- 
ply the  tangential  fields  — hence  the  name  “edge  elements”. 
For  higher  order  elements,  additional  degrees  of  freedom  on 
faces  must  be  introduced,  and  the  name  “vector  elements” 


has  now  largely  supplanted  “edge  elements”.  For  CT/LN  el- 
ements, some  researchers  have  associated  the  degree  of  free- 
dom with  the  tangential  field  at  the  centre  of  each  edge  [11]. 

B.  Vector  vs.  mixed-order  elements 

It  is  not  always  appreciated  that  being  of  mixed-order  is 
not  an  essential  property  of  vector  elements  per  se.  Com- 
plete sets  of  vector  elements  have  also  been  described  [9], 
with  degrees  of  freedom  proportional  to  tangential  field  com- 
ponents, as  for  mixed-order  elements.  (This  permits  explicit 
enforcement  of  tangential  field  continuity  only,  as  for  mixed- 
order  elements.  As  is  well  known,  this  type  of  field  continu- 
ity is  very  difficult  to  arrange  in  general  with  nodal-based  el- 
ements, which  are  also  generally  complete).  However,  such 
complete  sets  of  vector  element  produce  “wasted”  d.o.f.’s  for 
wave  eigenvalue  problems.  See  [12]  for  a comprehensive 
discussion  of  this.  In  essence,  Nedelec’s  constraints  provide 
mixed-order  elements  that  model  the  curl-space  as  efficiently 
as  possible,  for  a given  number  of  degrees  of  freedom.  Re- 
cent work  by  Webb  [ 1 3]  has  indicated  that  some  vector  elec- 
tromagnetic problems  are  more  efficiently  analyzed  using 
complete-order  vector  elements,  typically  when  the  solution 
is  dominated  by  electric  fields  strongly  “gradient”  in  nature. 
(The  specific  example  Webb  uses  is  an  iris  in  a waveguide, 
where  the  solution  is  strongly  dominated  by  quasi-static  field 
components). 

C.  Higher  order  elements 

Although  extending  the  “edge”  elements  to  higher  order 
became  a topic  of  interest  as  soon  as  the  CT/LN  elements 
achieved  widespread  acceptance,  it  remains  a topic  of  active 
research  at  present,  a decade  or  more  later.  Development  of 
such  elements  raises  a number  of  issues,  including:  hierar- 
chal  vs.  interpolator  behaviour;  methods  for  the  construc- 
tion of  the  element  shape  functions;  the  interpretation  of  the 
degrees  of  freedom;  the  construction  of  prototype  elemental 
matrices  (analytical  vs.  quadrature);  and  the  efficient  itera- 
tive solution  of  the  poorly  conditioned  linear  algebra  systems 
which  unfortunately  often  result. 

C.  1 Hierarchal  higher  order  LT/QN  elements 

For  mesh  refinement/enrichment  purposes,  hierarchal  el- 
ements are  required,  and  this  paper  considers  only  the  use 
of  such  elements.  Interpolator  elements  have  been  com- 
prehensively described  in  [14].  Two  specific  hierarchal  ele- 
ments have  been  used;  the  work  was  originally  undertaken 
[15],  [16],  [7]  using  those  proposed  by  Savage  [17];  see  Ta- 
ble I.  Subsequently,  the  elements  proposed  by  Andersen  and 
Volakis  [18],  [19]  have  also  been  implemented. 

As  per  Nedelec’s  definitions  of  suitable  mixed  order  ele- 
ment, there  are  twenty  vector  based  functions  (v.b.f.’s)  and 
degrees  of  freedom  (d.o.f.’s)  per  tetrahedron.  The  lowest  or- 
der v.b.f.  — the  Whitney  form  — has  the  usual  properties, 
and  the  accompanying  d.o.f.  is  proportional  to  the  tangen- 
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CT/LN 

Edge-based 

1 per  edge 

LT/QN  | 

Edge-based 

1 per  edge 

V(  XiXj) 

Face-based 

2 per  face 
(out  of  3 possible) 

MAjVAfc  - AfeVAj) 

TABLE  I 

Savage’s  LT/QN  hierarchal  elements. 


tial  electric  field  on  the  edge.  Since  the  element  shape  func- 
tions are  hierarchal,  higher-order  d.o.f.’s  are  not  located  at 
specific  points,  but  defined  as  weighted  field  quantities  in- 
tegrated over  the  relevant  edge  or  face.  The  latter  has  the 
interpretation  as  the  flux  through  the  face. 

Many  other  hierarchal  elements  have  been  published,  in 
particular  of  LT/QN  order.  Most  of  these  (including  those  of 
Savage  described  above)  can  be  seen  as  variants  of  the  el- 
ements proposed  by  Webb  and  Forghani  [20].  (Indeed,  not 
only  are  these  variants  on  a theme,  they  are  also  linear  trans- 
forms, as  will  be  discussed  subsequently).  A number  are 
summarized  in  Table  II.  Note  that  all  the  face  elements  ex- 
clude (arbitrarily)  one  possible  combination  of  {i;j;  k};  this 
asymmetry  has  long  been  noted,  and  is  required  to  avoid  lin- 
early dependent  basis  functions. 

An  apparent  exception  to  this  are  the  elements  proposed 
by  Andersen  and  Volakis  [18],  [19];  the  additional  six  vec- 
tor based  functions  for  the  edges  are  apparently  of  quadratic 
order.  However,  the  Andersen  and  Volakis  elements  are  lin- 
ear transforms  of  the  (non-hierarchal)  elements  proposed  in 
[2 1 ] (the  explicit  transform  for  the  two  dimensional  case  was 
given  in  [19])  and  the  hierarchal  elements  proposed  by  Sav- 
age [17]  — and  used  here  — are  in  turn  linear  combinations 
of  those  in  [21],  thus  the  Andersen  and  Volakis  elements  are 
linear  transforms  of  Savage’s  [17]. 

It  might  seem  strange  that  the  Andersen  and  Volakis  el- 
ements, with  apparently  quadratic  behaviour,  can  be  ex- 
pressed as  linear  combinations  of  elements  with  at  most  lin- 
ear field  dependence.  This  is  a consequence  of  the  mixed- 
order  nature  of  the  basis  functions,  and  of  course  the  linearly 
dependent  nature  of  simplex  coordinates  (]T)ili  ^ i = l,with 
N = 2,  3 or  4 in  one,  two  or  three  dimensions  respectively). 
The  (AiVAj  — AjVAi)  term  is  of  course  the  Whitney  ele- 
ment, with  CT/LN  behaviour  along  edges;  multiplication  by 
the  (A i — A j)  term  (actually  the  Legendre  polynomial  Pi  re- 
defined on  the  interval  [0, 1])  yields  the  LT/QN  behaviour. 

These  elements  are  generally  constructed  by  “inspection”, 
using  the  properties  of  simplex  coordinates,  and  the  gradi- 
ents thereof.  Nedelec  required  that  these  functions  be  uni- 
solvent (that  is,  linearly  independent)  and  conforming  (that 
is,  d.o.f.’s  are  proportional  to  integrals  of  the  tangential  field 
along  edges,  or  over  faces).  The  latter  is  easily  established  by 
inspection,  but  the  former  is  less  obvious.  The  present  author 


CT/LN  - all 

| Edge-based 

1 per  edge 

AjVA,-  — XjVXi  | 

LT/QN  — Savage  | 

1 per  edge 

V(XiXj) 

2 per  face 

Ai(A,-VA*  -AfcVA.  ) 

(and  {J;  i;k]  but  not  {k\i;j}) 

HMH 

LT/QN  — Webb  and  Forghani 

Edge-based 

1 per  edge 

Face-based 

2 per  face 

AAfcVA; 

(and  {j;  k ; i } but  not  {i;j;k}) 

LT/QN  — Andersen  and  Volakis 

Edge-based 

1 per  edge 

(Ai  - Aj)x 
(A;VAj  — AjVAi) 

Face-based 

2 per  face 

(as  for  Savage’s  elements) 

- \kV\j) 

TABLE  II 

Comparison  of  various  hierarchal  LT/QN  element  schemes. 


has  shown  that  all  these  basis  functions  satisfy  the  Nedelec 
constraints  (restrictions  on  the  properties  of  the  polynomial 
spaces  from  which  they  are  chosen)  and  thus  (from  [8,  The- 
orem 1])  the  elements  are  indeed  both  conforming  and  uni- 
solvent. This  proof  requires  expressing  the  basis  functions 
in  Cartesian  coordinate  form  and  then  testing  the  Nedelec 
constraints  explicitly;  it  will  not  be  detailed  here. 

There  is  another  school  of  thought  regarding  the  con- 
struction of  higher-order  basis  functions,  which  might  be 
described  as  the  degree  of  freedom-centered  approach  (as 
opposed  to  the  above,  which  could  be  described  as  the  ba- 
sis function-centered  approach).  Salazar-Palma  et  al.  [11] 
use  elements  from  the  Nedelec  polynomial  space  and  en- 
force Lagrangian  interpolator  properties  on  the  degrees  of 
freedom.  This  produces  interpolator  elements  with  well- 
defined  degrees  of  freedom  at  points , but  this  is  not  possible 
in  general  with  higher-order  hierarchal  elements.  Yioultsis 
and  Tsikboukis  take  a similar  degree-of-freedom  centered 
approach,  but  working  with  simplex  instead  of  Cartesian  co- 
ordinates [22]. 

IV.  IMPLEMENTATION  ISSUES 
A . Finite  element  discretization 

The  finite  element  discretization  of  the  volumetric  integral 
term  is  identical  to  that  of  cavity  eigenanalysis.  This  has 
been  described  in  several  references  (such  as  [3],  [23],  [21], 
[24])  and  will  not  be  discussed  further  here. 

Discretization  of  the  surface  integral  terms,  which  arises 
due  to  the  introduction  of  the  ports,  requires  compatible  sur- 
face basis  functions.  This  is  discussed  in  Jin  in  detail  [3, 
§8.5],  and  need  only  be  outlined  here,  since  the  extension  to 
higher-order  elements  is  obvious.  Generation  of  the  volu- 
metric tetrahedral  mesh  automatically  generates  a triangular 
surface  mesh.  Suitable  basis  functions  also  implicitly  de- 
fined, as  follows: 


DAVIDSON:  HIGHER-ORDER  (LT/QN)  VECTOR  FINITE  ELEMENTS  FOR  WAVEGUIDE  ANALYSIS 


5 


fix  Es  — ^ §■  = {Es}t{Ss}  (10) 

i=  1 

with  { Es } are  the  degrees  of  freedom  associated  with  the 
(surface)  triangular  element. 

The  surface  basis  functions  are: 

5f  = nxNf  (11) 

Arf  are  the  appropriate  tetrahedral  functions  for  the  face. 

Note  that  for  CT/LN  elements,  this  magnetic  surface  cur- 
rent (n  x Es)  discretization  is  identical  to  that  produced  by 
standard  moment  method  RWG  elements  [25],  providing  a 
linear  tangential/constant  normal  approximation  of  the  cur- 
rent. 

B.  Elemental  matrices  and  matrix  assembly 

A new  elemental  matrix  and  vector  are  required: 

[Bs]  = J J 1Ss-SsdS  (12) 

{6s}  = IJ  Ss  ■ (Einc  x n)dS  (13) 

Here,  7 is  as  in  eqn.  (3),  and  Elnc  is  the  incident  field,  as 
before. 

[Bs]  can  be  evaluated  in  closed  form,  since  it  involves  the 
integral  of  simplex  coordinates  over  both  ports  — these  inte- 
grals are  known  analytically.  {bs}  requires  quadrature,  since 
it  involves  the  product  of  the  incident  mode,  typically  a si- 
nusoid or  product  of  sinusoids,  with  a vector  based  function. 
A four-point  symmetric  rule  [26]  was  generally  found  to  be 
sufficient,  although  a six-point  rule  was  also  implemented. 

The  system  matrix  [A]  is  assembled  from  [5],  [ T } and  [B]; 
the  forcing  vector  is  {bs},  resulting  in  the  conventional  lin- 
ear system  [A]{x}  = {bs}  with  {x}  the  vector  of  degrees 
of  freedom  to  be  solved  for.  All  these  terms  are  frequency 
dependent,  and  [B]  and  {bs}  are  additionally  also  dependent 
on  the  mode  number  and/or  mode  type  (TE  or  TM),  either 
via  the  propagation  constant  or  the  modal  eigenfunction. 

Once  the  system  matrix  and  right  hand  side  vector  have 
been  assembled,  the  system  is  solved  (for  multiple  RHS’s, 
if  the  full  S-matrix  is  required)  and  the  S parameters  are  ex- 
tracted as  already  discussed  in  Section  II-B. 

The  [S]  and  [T]  elemental  matrices  may  be  pre-computed, 
using  explicit  forms  as  given  in  [3],  [23],  [21],  [24],  [15]. 
However,  a non-trivial  amount  of  analytical  work  is  re- 
quired for  new  elements,  and  the  use  of  cubature  (three- 
dimensional  quadrature)  permits  far  quicker  program  devel- 
opment; new  element  basis  functions  (and  their  gradients, 
which  are  straightforward  to  compute  analytically),  can  be 
added  in  very  quickly.  Since  the  functions  being  integrated 
are  polynomials,  and  very  efficient  rules  exist  for  integration 
of  polynomials  over  simplexes,  the  computational  overhead 


h 

(mm) 

N 

CT/LN 

LT/QN 

D.o.f.'s 

ls2ll 

(dB> 

LS2l 

(°> 

D.o.f.'s 

15-211 

(dB) 

LS21 

(°> 

20.4 

20 

14 

• 16.18 

134.4 

96 

-0,4113 

22.11 

13.8 

43 

29 

-4.538 

5.602 

202 

-0.3001 

6.077 

9.87 

I0S 

73 

-1.157 

-38.57 

498 

-0.01575 

-0.2293 

7.31 

234 

172 

-0.5191 

-19.72 

1 144 

-2.78E-3 

0.0448 

6.47 

369 

273 

-0.7414 

-31.87 

1810 

-2.26E-3 

-0,1619 

4.99 

697 

600 

-0.1944 

-3.615 

3700 

-4.34E-4 

0.0251 

4.22 

1123 

962 

-0.1455 

-0.561 

5948 

-9.56E-4 

0.0546 

TABLE  III 

521  FOR  an  empty  section  of  waveguide  for  the  dominant 
TE10  mode  as  a function  of  the  average  edge  length  h [mm] 
and  number  of  elements  N. 


is  modest,  and  of  O(N).  As  an  example,  a symmetric  rule 
of  degree  of  precision  four  requires  only  eleven  points,  and 
this  is  sufficient  for  LT/QN  basis  functions.  (The  elemental 
matrix  entries  in  this  case  require  at  most  the  integration  of 
the  product  of  quadratics,  i.e.  a polynomial  of  order  four). 

V.  Results 

The  theory  discussed  above  has  been  implemented  in  a 
finite  element  code  developed  by  the  author,  his  students  and 
industrial  colleagues.  The  code  uses  the  same  graphical  user 
interface  as  the  commercial  package  FEKO,  and  is  called 
FEMFEKO  [27]. 

A.  Empty  guide 

An  empty  section  of  waveguide  provides  a useful  test  of 
the  performance  of  the  elements,  since  the  results  are  known 
exactly.  In  Table  III,  the  transmission  coefficient  of  a hollow 
piece  of  X-band  waveguide,  40mm  long,  is  presented.  The 
analytical  solution  is  trivial;  the  transmission  coefficient  is  1 
(OdB),  with  phase  angle  0°  in  Jin’s  formulation,  or  — kZl0£ 
for  the  extended  formulation  presented  in  this  paper,  as  dis- 
cussed in  Section  II.  (The  results  in  Table  III  were  gener- 
ated with  Jin’s  original  formulation,  hence  the  zero  phase 
angles).  The  20  element  result  is  of  course  highly  inaccurate, 
since  the  problem  is  badly  under-discretized  with  so  few  el- 
ements. (The  guide  wavelength  was  48.630mm).  The  mesh 
refinement  used  in  Table  III  was  a simple  h-uniform  scheme. 

Eigenvalue  problems  are  appealing  since  one  quantity  (the 
eigenvalue)  can  be  checked  for  convergence  (and  it  is  also 
known  that  the  eigenvalue  is  variational  [28]);  for  example, 
Savage  and  Peterson  reported  convergence  results  for  LT/QN 
elements  for  eigenvalue  problems  in  [21].  Investigating  the 
convergence  of  scattering  parameters  is  somewhat  more  dif- 
ficult. The  energy  conservation  term,  |Sn|2  + IS21J2,  is  a 
useful  overall  solution  quality  indicator,  suggested  by  Jin, 
and  will  be  used  here.  For  a lossless  structure,  this  should 
of  course  be  unity.  Results  are  presented  in  Figs.l  and  2. 
The  former  shows  a consistently  lower  error  (result  closer 
to  unity)  for  the  same  h — i.e.  the  same  mesh  — and  thus 
modelling  and  pre-processing  effort,  although  of  course  the 
solution  using  LT/QN  elements  uses  many  more  degrees  of 
freedom.  The  latter  shows  a consistently  lower  error  for 
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Fig.  1.  |Sn|2  + 1 521 12  versus  h,  the  average  edge  length  in  the  mesh. 


Fig.  2.  |5n  |2  + |52i  I2  versus  the  number  of  degrees  of  freedom  used. 


the  same  number  of  degrees  of  freedom  — and  thus  com- 
putational effort  as  well,  presuming  that  efficient  solvers  are 
available  and  unaffected  by  the  use  of  the  higher-order  ele- 
ments (an  assumption  that  will  be  discussed  subsequently). 
This  indeed  is  the  major  motivation  for  using  higher-order 
elements. 

B.  E plane  bend 

As  a test  of  the  general  formulation,  and  also  of  the  rel- 
ative performance  of  the  elements,  an  E-plane  bend  will  be 
analyzed  (refer  to  Fig.  3).  The  analysis  will  be  performed 
for  X-band  (±  8-12  GHz)  waveguide.  It  will  be  seen  that  at 
the  lower  frequency  band,  the  bend  is  largely  transparent,  but 
towards  the  upper  end  of  the  frequency  band,  the  effects  of 
the  bend  become  significant.  This  problem  has  an  (approxi- 
mate) analytical  solution,  which  was  first  derived  some  sixty 
years  ago  by  Marcuvitz,  Schwinger  and  colleagues,  and  sub- 
sequently documented  in  [1]. 

This  problem  was  modelled  with  a section  of  “dummy” 
waveguide,  to  ensure  that  only  the  dominant  TEio  mode  is 
present  at  the  ports.  A half-length  of  40mm  is  sufficient; 


this  translates  into  around  30mm  of  waveguide  between  the 
bend  and  the  port  (that  is,  £ — b = 30mm).  The  geometry  is 
shown  in  Fig.  3.  It  is  assumed  that  the  waveguide  is  air- filled, 
approximated  as  free-space. 


Fig.  3.  The  E-plane  waveguide  bend.  Total  (half)  length  of  the  bend  is  L 


Fig.  4.  Side  view  of  the  E-plane  waveguide  bend. 


Fig.  5.  Equivalent  circuit  of  the  E-plane  waveguide  bend. 


B.  1 Comparing  with  Marcuvitz’s  approximate  analytical  re- 
sults 

The  equivalent  circuit  for  the  E-plane  bend  of  Fig.  3 as 
derived  by  Marcuvitz  [1,  pp.312ff]  is  pure  susceptance  — 
a shunt  inductor,  — jB;  (B  > 0)  — at  terminal  planes  T, 
located  distance  d from  the  outer  comer  of  the  bend.  See 
Figs.  4 and  5.  In  the  following,  Yq(=  1 jZf)  is  the  waveguide 
characteristic  admittance;  Z0  = 77/ 1 — (A/(2a))2  with 
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A the  free-space  wavelength,  r/  the  wave  impedance  of  free 
space  and  Xg  = X/y/l  - (A/(2 a))2  the  guide  wavelength. 

The  formulation  presented  in  this  paper  computes  5 pa- 
rameters, whereas  this  model  is  given  in  terms  of  a shunt  sus- 
ceptance.  There  are  several  methods  that  may  be  used  to  con- 
vert Marcuvitz’s  model  to  the  same  format  as  the  present  for- 
mulation. A straightforward  technique  is  to  find  the  equiv- 
alent ABCD  model  for  the  shunt  susceptance,  using  tables 
available  in  standard  texts  (for  example,  [29]);  then  convert 
this  to  5 parameters;  and  finally,  embed  this  within  a section 
of  guide  t — d.  This  procedure  will  now  be  applied. 

Firstly,  the  normalized  shunt  susceptance  and  distance  d 
can  be  read  off  [1,  Fig. 5. 28-3].  This  is  then  multiplied  by  the 
normalized  factor  2Y0b/Xg , giving  f3ShUnt- 

The  ABCD  matrix  [29,  Table  4.1]  for  this  shunt  load  is: 

A = 1;  B = 0;  C = -jBshunu  D=  1 (14) 

The  S-parameters  are  [29,  Table  4.2]: 

A + B/Z0  -CZo-D 
A + B/Z0  4-  CZo  4-  D 
2 (AD  - BC ) 

A + B/Z0  4-  CZq  T D 

2 

A + B j Zq  -f-  C Zo  4"  D 
-A  + B/Zo~CZo  + D 
A 4"  B / Zo  4"  CZo  4 ~ D 

Finally,  the  shunt  load  is  embedded  in  a line  of  length 
d — L This  amounts  to  changing  the  phase  of  the  indi- 
vidual S-parameters  by  e~^20 , with  0 = kZl0{£  - d)  [29, 
p.202-204],  and  kZlQ  = 2ir/Xg  the  wavenumber  of  the  TEio 
mode.  This  now  permits  direct  comparison  between  the  re- 
sults computed  using  the  FEM  and  the  approximate  results 
derived  by  Marcuvitz. 

B.2  Results 

For  the  results  to  be  presented,  two  meshes  were  gen- 
erated; a “coarse”mesh  with  an  average  edge  length  of 
about  6mm,  and  a “fine’'  mesh,  with  3.5mm  average  edge 
length.  The  problem  was  run  over  a frequency  range  of  8.25- 
12.25  GHz,  with  an  accompanying  guide  wavelength  vary- 
ing from  about  60-29mm.  Thus,  at  the  highest  frequency, 
the  coarse  mesh  was  about  A/5  — too  coarse  for  the  CT/LN 
elements  to  generate  reliable  solutions.  The  fine  mesh  should 
be  satisfactory. 

Results  are  presented  in  Figs.  6 and  7 for  the  CT/LN  ele- 
ments for  the  coarse  elements.  The  Sn  results  for  the  coarse 
mesh  are  indeed  very  inaccurate.  However,  the  fine  mesh 
yields  acceptable  results,  although  Sn  is  still  not  very  accu- 
rately computed.  (Note  that  at  the  low  frequency  end,  Sn  is 
small,  and  a very  accurate  solution  will  be  required  to  obtain 
good  agreement.) 


Clearly,  this  is  case  where  the  LT/QN  elements  will  be 
required  to  obtain  really  good  results.  These  are  presented  in 
Figs.  8 and  9 for  the  LT/QN  elements.  In  this  case,  the  phase 
comparison  is  also  shown,  in  Figs.  10  and  11.  Even  with 
the  coarse  mesh,  the  results  are  now  acceptable;  for  the  fine 
mesh,  the  agreement  is  excellent  across  the  entire  frequency 
band. 

Both  Savage’s  and  Andersen  and  Volakis’s  LT/QN  ele- 
ments were  used;  both  perform  very  well  from  a viewpoint  of 
accuracy.  There  was  no  discernible  difference  between  the 
S-parameter  results  computed  using  them.  Since  the  basis 
functions  are  related  by  a linear  transformation,  this  is  to  be 
expected.  However,  both  generated  ill-conditioned  matrices 
on  occasion.  This  remains  a problem  and  will  be  discussed 
in  the  next  section. 

C.  Iterative  solver  convergence  and  computational  effi- 
ciency 

A problem  which  has  not  been  addressed  in  this  paper 
is  the  issue  of  computational  efficiency . Unfortunately,  the 
higher-order  elements  appear  to  generate  ill-conditioned  ma- 
trices. The  results  presented  here  were  obtained  using  iter- 
ative solvers  for  the  linear  algebra  — variants  on  the  conju- 
gate gradient  scheme  (CG,  Bi-CG),  QMR  and  GMRES,  with 
simple  diagonal  pre-conditioning,  where  relevant  — but  all 
converged  erratically,  some  schemes  converging  rapidly  at 
certain  frequencies,  and  then  converging  very  slowly  at  oth- 
ers, and  also  exhibiting  different  convergence  behaviour  for 
different  geometries.  {All  converged  at  an  acceptable  pace 
for  the  CT/LN  elements).  This  is  a problem  which  is  only 
hinted  at  in  much  of  the  literature  on  higher-order  vector 
based  elements.  An  exception  is  the  work  by  Webb  [13]; 
he  proposes  a scheme  to  improve  the  matrix  conditioning 
by  at  least  partially  orthogonalizing  the  higher-order  basis 
functions.  Other  recent  approaches  have  focussed  on  the  use 
of  more  sophisticated  pre-conditioners.  Incomplete  LU  pre- 
conditioning is  one  possibility;  another  is  the  use  of  a di- 
rect solution  of  the  CT/LN  solution  (which  can  generally  be 
computed  quite  cheaply)  as  a pre-conditioner  for  the  LT/QN 
matrix. 

VI.  Conclusions 

This  paper  has  discussed  the  use  of  LT/QN  elements  for 
waveguide  analysis.  As  would  be  expected,  the  LT/QN  el- 
ements give  much  better  solutions  for  the  same  mesh  than 
CT/LN  elements;  this  remains  true  if  the  number  of  degrees 
of  freedom  are  compared.  Which  of  the  many  published 
LT/QN  elements  are  used  appears  insignificant  in  terms  of 
solution  accuracy  (at  least  for  the  E-plane  bend  analyzed  in 
this  paper,  as  well  as  several  other  problems  not  reported 
here);  the  choice  of  element  does  impact  on  the  convergence 
of  the  iterative  solver,  but  unfortunately  not  in  a consistent 
fashion. 

Although  the  performance  of  higher- order  elements  is 
usually  compared  with  that  of  lower-order  elements  in  terms 


5n  = 
S\2  ~ 

s2 1 = 

s22  = 
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Fig.  6.  Magnitude  of  the  reflection  coefficient  for  the  E-plane  bend  in  X 
band  guide;  CT/LN  elements. 


Fig.  9.  Magnitude  of  the  transmission  coefficient  for  the  E-plane  bend  in 
X-band  guide;  LT/QN  elements. 


Fig.  7.  Magnitude  of  the  transmission  coefficient  for  the  E-plane  bend  in  Fig.  10.  Phase  of  the  reflection  coefficient  for  the  E-plane  bend  in  X-band 
X-band  guide;  CT/LN  elements.  guide;  LT/QN  elements. 


Fig.  8.  Magnitude  of  the  reflection  coefficient  for  the  E-plane  bend  in  X- 
band  guide;  LT/QN  elements. 


Fig.  11.  Phase  of  the  transmission  coefficient  for  the  E-plane  bend  in  X- 
band  guide;  LT/QN  elements. 
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of  the  number  of  degrees  of  freedom  required  for  a partic- 
ular accuracy,  it  is  worth  making  the  point  that  the  geomet- 
rical pre-processing  (and  to  a lesser  extent,  post-processing) 
required  in  a real-world  FE  code  is  largely  a function  of  the 
number  of  elements,  rather  than  of  number  of  the  degrees  of 
freedom.  The  time  required  for  this  can  become  a significant 
fraction  of  the  total  run-time  of  the  code.  This  is  another 
practical  advantage  of  higher-order  elements  not  often  men- 
tioned in  the  literature. 

Yet  higher-order  schemes  — quadratic  tangential  /cubic 
normal  and  also  cubic  tangential  / quartic  normal  basis  func- 
tions have  been  published  (eg  [21],  [17]);  the  extension  of 
the  present  work  to  these  is  relatively  straightforward  the- 
oretically, but  implementing  these  will  require  meticulous 
work. 

From  the  viewpoint  of  waveguide  applications,  further 
waveguide  implementations  could  include  TM  modes;  this 
would  require  an  H solver  (due  to  the  boundary  condition) 
if  the  same  formulation  is  used,  but  this  should  present  no 
major  problems.  Multi-mode  analysis  is  already  included  in 
the  formulation  (and  been  implemented  in  the  code  which 
generated  these  results),  but  had  not  been  tested  at  the  time 
of  writing. 

Extending  the  formulation  to  right-circular  cylindrical 
waveguides  (or  indeed  any  shape  for  which  the  eigenmodes 
are  known  in  closed  form  and  are  transverse  in  nature)  would 
be  moderately  straightforward.  Extensions  to  include  arbi- 
trary, numerically-determined  modes  would  be  a desirable 
future  addition.  Inhomogeneously  loaded  waveguide  poses 
new  challenges,  since  one  needs  to  simultaneously  solve  for 
the  tangential  and  axial  fields.  Hybrid  schemes,  using  vec- 
tor elements  for  the  former  and  nodal  schemes  for  the  latter 
have  been  used  successfully  [30],  but  there  are  a variety  of 
approaches  to  this  requiring  exploration.  These  also  pose 
some  problems  for  the  present  formulation,  which  relies  on 
the  transverse  nature  of  the  eigenmodes  to  obtain  the  bound- 
ary condition  at  the  ports. 

Another  possible  extension  would  be  to  include  anisotropic 
materials;  the  formulation  for  diagonally  anisotropic  materi- 
als for  CT/LN  elements  is  available  [24]  and  this  would  have 
to  be  extended  to  fully  anisotropic  media  and  LT/QN  ele- 
ments. Non-linear  materials  are  however  probably  best  left 
to  a time-domain  approach,  in  particular  the  FDTD. 
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